%% RCT2Scale empirical analysis
% Replicates Figure 10 left & center columns, Figure A-11, Table A-4, Figure A-12

%% Parametrization
% Number of bins for PIT histograms
bin_num = 5;

% Set max iteration number and tolerance
max_iter = 5000;
tol = 10^-6;

% Random seed
seed = 101;
rng(seed)

% Sample size of simulated coverage & S stats
sim_cov = 5000;
sim_S = 5000;

%% Case 1: nu=0 (grouping 1)
% Initialization
N = 66;
J = 2;

% Extract study estimates and ses
sheet1 = 'Estimates';
sheet2 = 'SE';
sheet3 = 'Sample Size';
estimates1 = xlsread('../Data/RCT2Scale_data.xlsx', sheet1, 'A2:C67'); 
se1 = xlsread('../Data/RCT2Scale_data.xlsx', sheet2, 'A2:C67');   
size1 = xlsread('../Data/RCT2Scale_Data.xlsx', sheet3, 'A2:C67'); 

% Randomly assign baseline and validation studies
[hat_theta_1, hat_vartheta_1, rep_base_se2_1, rep_val_se2_1, size_base_1, size_val_1] = random_assignment(estimates1, se1, size1, seed);

% Compute posterior means and variances for tau
[bar_tau_1, bar_V_tau_1] = compute_posterior_moments(hat_theta_1, rep_base_se2_1, 0);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[predict_se_1, empirical_cov_1, cov_p_value_1] = check_CI(hat_vartheta_1, rep_val_se2_1, bar_tau_1, bar_V_tau_1, 0, cov_stat_sim);
xtick = 0:10:70;
ytick = -25:5:30;
[interval_1, order_1] = plot_intervals(hat_vartheta_1, bar_tau_1, predict_se_1, xtick, ytick);
saveas(interval_1, "../Results/Benchmark_Interval_N_66_J_2_g1.png")

% Generate the S statistics and p-value
S_stat_sim_1 = simulate_S_stat_distribution(N, sim_S);
ytick = 0:0.1:0.5;
[S_stat_1, S_p_value_1, PIT_plot_1] = PITs(hat_vartheta_1, bar_tau_1, predict_se_1, S_stat_sim_1, bin_num, ytick);
saveas(PIT_plot_1, "../Results/Benchmark_PIT_N_66_J_2_g1.png")

%% Case 2: EB (grouping 1)

% Estimate nu from pooled data
EB_est_nu_2 = MLE_nu(hat_theta_1, rep_base_se2_1, 0, max_iter, tol);

% Compute posterior means and variances for tau
[bar_tau_2, bar_V_tau_2] = compute_posterior_moments(hat_theta_1, rep_base_se2_1, EB_est_nu_2);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[predict_se_2, empirical_cov_2, cov_p_value_2] = check_CI(hat_vartheta_1, rep_val_se2_1, bar_tau_2, bar_V_tau_2, EB_est_nu_2, cov_stat_sim);
xtick = 0:10:70;
ytick = -25:5:30;
interval_2 = plot_intervals_fixord(hat_vartheta_1, bar_tau_2, predict_se_2, order_1, xtick, ytick);
saveas(interval_2, "../Results/EB_Interval_N_66_J_2_g1.png")

% Generate the S statistics and p-value
S_stat_sim_2 = simulate_S_stat_distribution(N, sim_S);
ytick = 0:0.1:0.6;
[S_stat_2, S_p_value_2, PIT_plot_2] = PITs(hat_vartheta_1, bar_tau_2, predict_se_2, S_stat_sim_2, bin_num, ytick);
saveas(PIT_plot_2, "../Results/EB_PIT_N_66_J_2_g1.png")

% Plot the ratio of nu to se
yscale_rat = [0, 120];
yscale_rat_asym = [0, 0.7];
nu_se_fig_2 = plot_nu_se_ratio(EB_est_nu_2, rep_base_se2_1, rep_val_se2_1, yscale_rat, xtick);
nu_asymse_fig_2 = plot_nu_asymse_ratio(EB_est_nu_2, rep_base_se2_1, rep_val_se2_1, size_base_1, size_val_1, yscale_rat_asym, xtick);
saveas(nu_se_fig_2, "../Results/nu_se_ratio_N_66_J_2_g1.png")
saveas(nu_asymse_fig_2, "../Results/nu_asymse_ratio_N_66_J_2_g1.png")

%% Case 3: nu=0 (grouping 2)
% Initialization
N = 66;
J = 2;

% Extract study estimates and ses
sheet1 = 'Estimates';
sheet2 = 'SE';
estimates3 = xlsread('../Data/RCT2Scale_data2.xlsx', sheet1, 'A2:C67'); 
se3 = xlsread('../Data/RCT2Scale_data2.xlsx', sheet2, 'A2:C67');   
size3 = xlsread('../Data/RCT2Scale_data2.xlsx', sheet3, 'A2:C67'); 

% Randomly assign baseline and validation studies
[hat_theta_3, hat_vartheta_3, rep_base_se2_3, rep_val_se2_3, size_base_3, size_val_3] = random_assignment(estimates3, se3, size3, seed);

% Compute posterior means and variances for tau
[bar_tau_3, bar_V_tau_3] = compute_posterior_moments(hat_theta_3, rep_base_se2_3, 0);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[predict_se_3, empirical_cov_3, cov_p_value_3] = check_CI(hat_vartheta_3, rep_val_se2_3, bar_tau_3, bar_V_tau_3, 0, cov_stat_sim);
xtick = 0:10:70;
ytick = -25:5:20;
[interval_3, order_3] = plot_intervals(hat_vartheta_3, bar_tau_3, predict_se_3, xtick, ytick);
saveas(interval_3, "../Results/Benchmark_Interval_N_66_J_2_g2.png")

% Generate the S statistics and p-value
S_stat_sim_3 = simulate_S_stat_distribution(N, sim_S);
ytick = 0:0.1:0.5;
[S_stat_3, S_p_value_3, PIT_plot_3] = PITs(hat_vartheta_3, bar_tau_3, predict_se_3, S_stat_sim_3, bin_num, ytick);
saveas(PIT_plot_3, "../Results/Benchmark_PIT_N_66_J_2_g2.png")

%% Case 4: EB (grouping 2)
% Estimate nu from pooled data
EB_est_nu_4 = MLE_nu(hat_theta_3, rep_base_se2_3, 0, max_iter, tol);

% Compute posterior means and variances for tau
[bar_tau_4, bar_V_tau_4] = compute_posterior_moments(hat_theta_3, rep_base_se2_3, EB_est_nu_4);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[predict_se_4, empirical_cov_4, cov_p_value_4] = check_CI(hat_vartheta_3, rep_val_se2_3, bar_tau_4, bar_V_tau_4, EB_est_nu_4, cov_stat_sim);
xtick = 0:10:70;
ytick = -25:5:20;
interval_4 = plot_intervals_fixord(hat_vartheta_3, bar_tau_4, predict_se_4, order_3, xtick, ytick);
saveas(interval_4, "../Results/EB_Interval_N_66_J_2_g2.png")

% Generate the S statistics and p-value
S_stat_sim_4 = simulate_S_stat_distribution(N, sim_S);
ytick = 0:0.1:0.6;
[S_stat_4, S_p_value_4, PIT_plot_4] = PITs(hat_vartheta_3, bar_tau_4, predict_se_4, S_stat_sim_4, bin_num, ytick);
saveas(PIT_plot_4, "../Results/EB_PIT_N_66_J_2_g2.png")

% Plot the ratio of nu to se
yscale_rat = [0, 120];
yscale_rat_asym = [0, 0.7];
nu_se_fig_4 = plot_nu_se_ratio(EB_est_nu_4, rep_base_se2_3, rep_val_se2_3, yscale_rat, xtick);
nu_asymse_fig_4 = plot_nu_asymse_ratio(EB_est_nu_4, rep_base_se2_3, rep_val_se2_3, size_base_3, size_val_3, yscale_rat_asym, xtick);
saveas(nu_se_fig_4, "../Results/nu_se_ratio_N_66_J_2_g2.png")
saveas(nu_asymse_fig_4, "../Results/nu_asymse_ratio_N_66_J_2_g2.png")

%% Summarize the benchmark results
empirical_cov_list = [empirical_cov_1, empirical_cov_2, empirical_cov_3, empirical_cov_4]
cov_p_value_list = [cov_p_value_1, cov_p_value_2, cov_p_value_3, cov_p_value_4]
S_stat_list = [S_stat_1, S_stat_2, S_stat_3, S_stat_4]
S_p_value_list = [S_p_value_1, S_p_value_2, S_p_value_3, S_p_value_4]